remove(list = ls())
##Install Development Version for Lavaan and semTools
##library("devtools")
##install.packages("lavaan", repos="http://www.da.ugent.be", type="source")
##devtools::install_github("simsem/semTools/semTools")

library(lavaan)
library(lavaan.survey)
library(OrthoPanels)
library(semTools)
library(foreign)
library(dplyr)
library(readstata13)

location <- "/Users/dillonlaaker/Box Sync/Stability/"
setwd(paste(location, "data/datasets", sep = ""))
liss <- read.dta13("liss_data.dta")
bes <- read.dta("bes_data.dta")
ncp <- read.dta("ncp_data.dta")

##LISS

model1 <- "imm8 =~ imm5_8 + imm2_8 + imm3_8 + imm4_8 + imm1_8 + imm6_8
imm9 =~ imm5_9 + imm2_9 + imm3_9 + imm4_9 + imm1_9 + imm6_9
imm10 =~ imm5_10 + imm2_10 + imm3_10 + imm4_10 + imm1_10 + imm6_10
imm11 =~ imm5_11 + imm2_11 + imm3_11 + imm4_11 + imm1_11 + imm6_11
imm12 =~ imm5_12 + imm2_12 + imm3_12 + imm4_12 + imm1_12 + imm6_12
imm13 =~ imm5_13 + imm2_13 + imm3_13 + imm4_13 + imm1_13 + imm6_13
imm14 =~ imm5_14 + imm2_14 + imm3_14 + imm4_14 + imm1_14 + imm6_14
imm16 =~ imm5_16 + imm2_16 + imm3_16 + imm4_16 + imm1_16 + imm6_16
imm17 =~ imm5_17 + imm2_17 + imm3_17 + imm4_17 + imm1_17 + imm6_17
imm9 ~ imm8
imm10 ~ imm9
imm11 ~ imm10
imm12 ~ imm11
imm13 ~ imm12
imm14 ~ imm13
imm16 ~ imm14
imm17 ~ imm16"
fit1 <- sem(model1, se="robust.huber.white", test="Satorra-Bentler", data=liss)
reliabl1 <- reliability(fit1)
est <- fitMeasures(fit1, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model1results <- parameterEstimates(fit1)
model1results <- model1results[c(55,56,57,58,59,60,61,62),c(4,5)]
model1results <- model1results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 8)
par2 <- rep(")", 8)
model1results$se <- paste(par1, model1results$se, par2, sep="")
model1results$est <- paste(model1results$est, model1results$se)
model1results$se <- NULL
coefnames <- c("$\\beta_{1,2}$", "$\\beta_{2,3}$", "$\\beta_{3,4}$", "$\\beta_{4,5}$", "$\\beta_{5,6}$", "$\\beta_{6,7}$", "$\\beta_{7,8}$", "$\\beta_{8,9}$")
model1results <- data.frame(coefnames, model1results)
coefnames1 <- c("$\\chi^{2}$", "df", "p-value", "CFI", "TLI", "RMSEA", "SRMSR", "AIC")
fitmeasures <- data.frame(est)
est <- fitmeasures %>% mutate(est= format(round(est, 2), nsmall = 2))
fitmeasures <- data.frame(coefnames1, est)
colnames(fitmeasures) <- c("coefnames", "est")
obs <- nobs(fit1)
obs <- data.frame("N", obs)
colnames(obs) <- c("coefnames", "est")
results1 <- rbind(model1results, fitmeasures, obs)

model2 <- "imm8 =~ imm5_8 + imm2_8 + imm3_8 + imm4_8 + imm1_8 + imm6_8
imm9 =~ imm5_9 + imm2_9 + imm3_9 + imm4_9 + imm1_9 + imm6_9
imm10 =~ imm5_10 + imm2_10 + imm3_10 + imm4_10 + imm1_10 + imm6_10
imm11 =~ imm5_11 + imm2_11 + imm3_11 + imm4_11 + imm1_11 + imm6_11
imm12 =~ imm5_12 + imm2_12 + imm3_12 + imm4_12 + imm1_12 + imm6_12
imm13 =~ imm5_13 + imm2_13 + imm3_13 + imm4_13 + imm1_13 + imm6_13
imm14 =~ imm5_14 + imm2_14 + imm3_14 + imm4_14 + imm1_14 + imm6_14
imm16 =~ imm5_16 + imm2_16 + imm3_16 + imm4_16 + imm1_16 + imm6_16
imm17 =~ imm5_17 + imm2_17 + imm3_17 + imm4_17 + imm1_17 + imm6_17
imm1_8 ~~ imm1_9 + imm1_10 + imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_9 ~~ imm1_10 + imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_10 ~~ imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_11 ~~ imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_12 ~~ imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_13 ~~ imm1_14 + imm1_16 + imm1_17
imm1_14 ~~ imm1_16 + imm1_17
imm1_16 ~~ imm1_17
imm2_8 ~~ imm2_9 + imm2_10 + imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_9 ~~ imm2_10 + imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_10 ~~ imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_11 ~~ imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_12 ~~ imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_13 ~~ imm2_14 + imm2_16 + imm2_17
imm2_14 ~~ imm2_16 + imm2_17
imm2_16 ~~ imm2_17
imm3_8 ~~ imm3_9 + imm3_10 + imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_9 ~~ imm3_10 + imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_10 ~~ imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_11 ~~ imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_12 ~~ imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_13 ~~ imm3_14 + imm3_16 + imm3_17
imm3_14 ~~ imm3_16 + imm3_17
imm3_16 ~~ imm3_17
imm4_8 ~~ imm4_9 + imm4_10 + imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_9 ~~ imm4_10 + imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_10 ~~ imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_11 ~~ imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_12 ~~ imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_13 ~~ imm4_14 + imm4_16 + imm4_17
imm4_14 ~~ imm4_16 + imm4_17
imm4_16 ~~ imm4_17
imm5_8 ~~ imm5_9 + imm5_10 + imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_9 ~~ imm5_10 + imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_10 ~~ imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_11 ~~ imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_12 ~~ imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_13 ~~ imm5_14 + imm5_16 + imm5_17
imm5_14 ~~ imm5_16 + imm5_17
imm5_16 ~~ imm5_17
imm6_8 ~~ imm6_9 + imm6_10 + imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_9 ~~ imm6_10 + imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_10 ~~ imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_11 ~~ imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_12 ~~ imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_13 ~~ imm6_14 + imm6_16 + imm6_17
imm6_14 ~~ imm6_16 + imm6_17
imm6_16 ~~ imm6_17
imm2_8 ~~ imm4_8
imm5_8 ~~ imm6_8
imm4_8 ~~ imm6_8
imm2_8 ~~ imm3_8
imm2_8 ~~ imm6_8
imm1_9 ~~ imm2_9 + imm3_9 + imm4_9 + imm5_9 + imm6_9
imm2_9 ~~ imm3_9 + imm4_9 + imm5_9 + imm6_9
imm3_9 ~~ imm4_9 + imm5_9 + imm6_9
imm4_9 ~~ imm5_9 + imm6_9
imm5_9 ~~ imm6_9
imm1_10 ~~ imm2_10 + imm3_10 + imm4_10 + imm5_10 + imm6_10
imm2_10 ~~ imm3_10 + imm4_10 + imm5_10 + imm6_10
imm3_10 ~~ imm4_10 + imm5_10 + imm6_10
imm4_10 ~~ imm5_10 + imm6_10
imm5_10 ~~ imm6_10
imm1_11 ~~ imm2_11 + imm3_11 + imm4_11 + imm5_11 + imm6_11
imm2_11 ~~ imm3_11 + imm4_11 + imm5_11 + imm6_11
imm3_11 ~~ imm4_11 + imm5_11 + imm6_11
imm4_11 ~~ imm5_11 + imm6_11
imm5_11 ~~ imm6_11
imm1_12 ~~ imm2_12 + imm3_12 + imm4_12 + imm5_12 + imm6_12
imm2_12 ~~ imm3_12 + imm4_12 + imm5_12 + imm6_12
imm3_12 ~~ imm4_12 + imm5_12 + imm6_12
imm4_12 ~~ imm5_12 + imm6_12
imm5_12 ~~ imm6_12
imm1_13 ~~ imm2_13 + imm3_13 + imm4_13 + imm5_13 + imm6_13
imm2_13 ~~ imm3_13 + imm4_13 + imm5_13 + imm6_13
imm3_13 ~~ imm4_13 + imm5_13 + imm6_13
imm4_13 ~~ imm5_13 + imm6_13
imm5_13 ~~ imm6_13
imm1_14 ~~ imm2_14 + imm3_14 + imm4_14 + imm5_14 + imm6_14
imm2_14 ~~ imm3_14 + imm4_14 + imm5_14 + imm6_14
imm3_14 ~~ imm4_14 + imm5_14 + imm6_14
imm4_14 ~~ imm5_14 + imm6_14
imm5_14 ~~ imm6_14
imm1_16 ~~ imm2_16 + imm3_16 + imm4_16 + imm5_16 + imm6_16
imm2_16 ~~ imm3_16 + imm4_16 + imm5_16 + imm6_16
imm3_16 ~~ imm4_16 + imm5_16 + imm6_16
imm4_16 ~~ imm5_16 + imm6_16
imm5_16 ~~ imm6_16
imm5_17 ~~ imm6_17
imm2_17 ~~ imm3_17
imm2_17 ~~ imm6_17
imm2_17 ~~ imm4_17
imm3_17 ~~ imm6_17
imm5_17 ~~ imm4_17
imm9 ~ imm8
imm10 ~ imm9
imm11 ~ imm10
imm12 ~ imm11
imm13 ~ imm12
imm14 ~ imm13
imm16 ~ imm14
imm17 ~ imm16"

fit2<-sem(model2, se="robust.huber.white", test="Satorra-Bentler", data=liss)
reliabl2 <- reliability(fit2)
est <- fitMeasures(fit2, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model2results <- parameterEstimates(fit2)
model2results <- model2results[c(387,388,389,390,391,392,393,394),c(4,5)]
model2results <- model2results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 8)
par2 <- rep(")", 8)
model2results$se <- paste(par1, model2results$se, par2, sep="")
model2results$est <- paste(model2results$est, model2results$se)
model2results$se <- NULL
fitmeasures2 <- data.frame(est)
fitmeasures2 <- fitmeasures2 %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit2)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results2 <- rbind(model2results, fitmeasures2, obs)

lissresults <- cbind(results1, results2)
colnames(lissresults) <- c("", "Model (1)", "Model (2)")

##BES

model1 <- "imm1 =~ immcul1 + immwelfare1 + immecon1
imm2 =~ immcul2 + immwelfare2 + immecon2
imm3 =~ immcul3 + immwelfare3 + immecon3
imm4 =~ immcul4 + immwelfare4 + immecon4
imm7 =~ immcul7 + immwelfare7 + immecon7
imm8 =~ immcul8 + immwelfare8 + immecon8
imm10 =~ immcul10 + immwelfare10 + immecon10
imm11 =~ immcul11 + immwelfare11 + immecon11
imm2 ~ imm1
imm3 ~ imm2
imm4 ~ imm3
imm7 ~ imm4
imm8 ~ imm7
imm10 ~ imm8
imm11 ~ imm10"
fit1 <- sem(model1, se="robust.huber.white", test="Satorra-Bentler", data=bes)
reliabb1 <- reliability(fit1)
est <- fitMeasures(fit1, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model1results <- parameterEstimates(fit1)
model1results <- model1results[c(25,26,27,28,29,30,31),c(4,5)]
model1results <- model1results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 7)
par2 <- rep(")", 7)
model1results$se <- paste(par1, model1results$se, par2, sep="")
model1results$est <- paste(model1results$est, model1results$se)
model1results$se <- NULL
model1results <- data.frame(model1results)
fix <- data.frame("")
colnames(fix) <- c("est")
model1results <- rbind(model1results, fix)
fitmeasures <- data.frame(est)
fitmeasures <- fitmeasures %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit1)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results1 <- rbind(model1results, fitmeasures, obs)

model2 <- "imm1 =~ immcul1 + immwelfare1 + immecon1
imm2 =~ immcul2 + immwelfare2 + immecon2
imm3 =~ immcul3 + immwelfare3 + immecon3
imm4 =~ immcul4 + immwelfare4 + immecon4
imm7 =~ immcul7 + immwelfare7 + immecon7
imm8 =~ immcul8 + immwelfare8 + immecon8
imm10 =~ immcul10 + immwelfare10 + immecon10
imm11 =~ immcul11 + immwelfare11 + immecon11
immecon1 ~~ immecon2 + immecon3 + immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon2 ~~ immecon3 + immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon3 ~~ immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon4 ~~ immecon7 + immecon8 + immecon10 + immecon11
immecon7 ~~ immecon8 + immecon10 + immecon11
immecon8 ~~ immecon10 + immecon11
immecon10 ~~ immecon11
immcul1 ~~ immcul2 + immcul3 + immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul2 ~~ immcul3 + immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul3 ~~ immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul4 ~~ immcul7 + immcul8 + immcul10 + immcul11
immcul7 ~~ immcul8 + immcul10 + immcul11
immcul8 ~~ immcul10 + immcul11
immcul10 ~~ immcul11
immwelfare1 ~~ immwelfare2 + immwelfare3 + immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare2 ~~ immwelfare3 + immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare3 ~~ immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare4 ~~ immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare7 ~~ immwelfare8 + immwelfare10 + immwelfare11
immwelfare8 ~~ immwelfare10 + immwelfare11
immwelfare10 ~~ immwelfare11
immcul1 ~~ immwelfare1
immcul2 ~~ immwelfare2
immcul3 ~~ immwelfare3
immecon4 ~~ immcul4 + immwelfare4
immecon7 ~~ immcul7 + immwelfare7
immecon8 ~~ immcul8 + immwelfare8
immecon10 ~~ immcul10 + immwelfare10
immecon11 ~~ immcul11 + immwelfare11
imm2 ~ imm1
imm3 ~ imm2
imm4 ~ imm3
imm7 ~ imm4
imm8 ~ imm7
imm10 ~ imm8
imm11 ~ imm10"
fit2 <- sem(model2, se="robust.huber.white", test="Satorra-Bentler", data=bes)
reliabb2 <- reliability(fit2)
est <- fitMeasures(fit2, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model2results <- parameterEstimates(fit2)
model2results <- model2results[c(122,123,124,125,126,127,128),c(4,5)]
model2results <- model2results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 7)
par2 <- rep(")", 7)
model2results$se <- paste(par1, model2results$se, par2, sep="")
model2results$est <- paste(model2results$est, model2results$se)
model2results$se <- NULL
fix <- data.frame("")
colnames(fix) <- c("est")
model2results <- rbind(model2results, fix)
fitmeasures2 <- data.frame(est)
fitmeasures2 <- fitmeasures2 %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit2)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results2 <- rbind(model2results, fitmeasures2, obs)

besresults <- cbind(results1, results2)
colnames(besresults) <- c("Model (3)", "Model (4)")


##NCP
model1 <- "imm3 =~ imm2_3 + imm1_3
imm4 =~ imm2_4 + imm1_4
imm5 =~ imm2_5 + imm1_5
imm6 =~ imm2_6 + imm1_6
imm7 =~ imm2_7 + imm1_7
imm8 =~ imm2_8 + imm1_8
imm4 ~ imm3
imm5 ~ imm4
imm6 ~ imm5
imm7 ~ imm6
imm8 ~ imm7"
fit1 <- sem(model1, se="robust.huber.white", test="Satorra-Bentler", data=ncp)
reliabn1 <- reliability(fit1)
est <- fitMeasures(fit1, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model1results <- parameterEstimates(fit1)
model1results <- model1results[c(13,14,15,16,17),c(4,5)]
model1results <- model1results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 5)
par2 <- rep(")", 5)
model1results$se <- paste(par1, model1results$se, par2, sep="")
model1results$est <- paste(model1results$est, model1results$se)
model1results$se <- NULL
model1results <- data.frame(model1results)
m <- matrix("", ncol = 1, nrow = 3)
fix <- data.frame(m)
colnames(fix) <- c("est")
model1results <- rbind(model1results, fix)
fitmeasures <- data.frame(est)
est <- fitmeasures %>% mutate(est= format(round(est, 2), nsmall = 2))
fitmeasures <- data.frame(est)
obs <- nobs(fit1)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results1 <- rbind(model1results, fitmeasures, obs)

model2 <- "imm3 =~ imm2_3 + imm1_3
imm4 =~ imm2_4 + imm1_4
imm5 =~ imm2_5 + imm1_5
imm6 =~ imm2_6 + imm1_6
imm7 =~ imm2_7 + imm1_7
imm8 =~ imm2_8 + imm1_8
imm1_3 ~~ imm1_4 + imm1_5 + imm1_6 + imm1_7
imm1_4 ~~ imm1_5 + imm1_6 + imm1_7 + imm1_8
imm1_5 ~~ imm1_6 + imm1_7 + imm1_8
imm1_6 ~~ imm1_7 + imm1_8
imm1_7 ~~ imm1_8
imm2_3 ~~ imm2_4 + imm2_5 + imm2_6 + imm2_7 + imm2_8
imm2_4 ~~ imm2_5 + imm2_6 + imm2_7 + imm2_8
imm2_5 ~~ imm2_6 + imm2_7 + imm2_8
imm2_6 ~~ imm2_7 + imm2_8
imm2_7 ~~ imm2_8
imm4 ~ imm3
imm5 ~ imm4
imm6 ~ imm5
imm7 ~ imm6
imm8 ~ imm7"
fit2 <- sem(model2, se="robust.huber.white", test="Satorra-Bentler", data=ncp)
reliabn2 <- reliability(fit2)
est <- fitMeasures(fit2, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model2results <- parameterEstimates(fit2)
model2results <- model2results[c(42,43,44,45,46),c(4,5)]
model2results <- model2results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 5)
par2 <- rep(")", 5)
model2results$se <- paste(par1, model2results$se, par2, sep="")
model2results$est <- paste(model2results$est, model2results$se)
model2results$se <- NULL
m <- matrix("", ncol = 1, nrow = 3)
fix <- data.frame(m)
colnames(fix) <- c("est")
model2results <- rbind(model2results, fix)
fitmeasures2 <- data.frame(est)
fitmeasures2 <- fitmeasures2 %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit2)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results2 <- rbind(model2results, fitmeasures2, obs)

ncpresults <- cbind(results1, results2)
colnames(ncpresults) <- c("Model (7)", "Model (8)")

fix <- data.frame(c(rep("", times=17)))
table <- cbind(lissresults, fix, besresults, fix, ncpresults)
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
comment          <- list()
comment$pos      <- list()
comment$pos[[1]] <- 1
comment$pos[[2]] <- 2
comment$pos[[3]] <- 3
comment$pos[[4]] <- 4
comment$pos[[5]] <- 5
comment$pos[[6]] <- 6
comment$pos[[7]] <- 7
comment$pos[[8]] <- 8
comment$command  <- c(rep("[1.5ex] \n", times=7), "\\midrule \n")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE, add.to.row = comment,   hline.after = NULL, file=paste(location, "Draft/tables/table3a.tex", sep = ""))